; $Id: pOmA.pro,v 1.4 2010-09-02 11:53:56 pkapyla Exp $
;
;set_plot,'PS',/copy
;device,file="pOmA.eps"
;
if !d.name eq 'PS' then begin
  device,xsize=18,ysize=24,yoffset=3,/color,/encapsul,BITS=8
  !p.charthick=4 & !p.thick=4 & !x.thick=4 & !y.thick=4
end
;
restore,file='../cool128x256x128a1/stresses.sav'
;
nr=n_elements(r) & ntheta=n_elements(theta) & rsinth=fltarr(nr,ntheta)
;
for i=0,nr-1 do begin
  for j=0,ntheta-1 do begin
    rsinth(i,j)=r(i)*sin(theta(j))
  endfor
endfor
;
nlev=35 & frac=0.02
;
!p.multi=[0,3,2] & !p.charsize=2.7 & !x.margin=[3.75,0] & !y.margin=[3,0.5]
bar='!20!s!A$!n!r!6'
;
restore,file='../cool128x256x128a1/stresses.sav'
Om0=0.05 & Om=uzphi/rsinth+Om0
levm=max([abs(1.-min(Om/Om0)),abs(max(Om/Om0)-1)])
lev=1.-levm+findgen(nlev)/float(nlev-1)*2.*levm
contour,Om/Om0,x,y,/iso,nl=nlev,/fi,lev=lev
colorbar,range=[min(lev),max(lev)],pos=[0.15,0.7,0.17,0.85],/vert,ytickformat='(F5.2)',yticks=2,ytickv=[min(lev),1.,max(lev)],yaxis=0,char=2.5
;
uuh=fltarr(nr,ntheta,3) & uuh(*,*,0)=uxphi & uuh(*,*,1)=uyphi & uuh(*,*,2)=uzphi/x
dimension=size(uuh) & nx=dimension[1] & ny=dimension[2]
ux=findgen(nx,ny) & uy=findgen(nx,ny)
xcoord=findgen(nx,ny) & ycoord=findgen(nx,ny)
for ix=0,nx-1 do begin
  for iy=0,ny-1 do begin
    ux[ix,iy]=uuh[ix,iy,0]*sin(theta[iy])+uuh[ix,iy,1]*cos(theta[iy])
    uy[ix,iy]=uuh[ix,iy,0]*cos(theta[iy])-uuh[ix,iy,1]*sin(theta[iy])
    xcoord[ix,iy] = r[ix]*sin(theta[iy])
    ycoord[ix,iy] = r[ix]*cos(theta[iy])
  end
end
;
partvelvec,ux,uy,xcoord,ycoord,fraction=frac,/over,length=0.05,col=255
;
!p.charsize=1.4 & xyouts,0.6,0.83,'!6Run A1' & !p.charsize=2.7
!p.charsize=1.4 & xyouts,0.72,-0.83,''+bar+'!7X!6/!7X!6!d0!n!6' & !p.charsize=2.7
;
restore,file='../cool128x256x128a2/stresses.sav'
Om0=0.1 & Om=uzphi/rsinth+Om0 
levm=max([abs(1.-min(Om/Om0)),abs(max(Om/Om0)-1)])
lev=1.-levm+findgen(nlev)/float(nlev-1)*2.*levm
contour,Om/Om0,x,y,/iso,nl=nlev,/fi,lev=lev
colorbar,range=[min(lev),max(lev)],pos=[0.485,0.7,0.505,0.85],/vert,ytickformat='(F5.2)',yticks=2,ytickv=[min(lev),1.,max(lev)],yaxis=0,char=2.5
;
uuh=fltarr(nr,ntheta,3) & uuh(*,*,0)=uxphi & uuh(*,*,1)=uyphi & uuh(*,*,2)=uzphi/x
for ix=0,nx-1 do begin
  for iy=0,ny-1 do begin
    ux[ix,iy]=uuh[ix,iy,0]*sin(theta[iy])+uuh[ix,iy,1]*cos(theta[iy])
    uy[ix,iy]=uuh[ix,iy,0]*cos(theta[iy])-uuh[ix,iy,1]*sin(theta[iy])
  end
end
;
partvelvec,ux,uy,xcoord,ycoord,fraction=frac,/over,length=0.05,col=255
;
!p.charsize=1.4 & xyouts,0.6,0.83,'!6Run A2' & !p.charsize=2.7
;
restore,file='../cool128x256x128a3/stresses.sav'
Om0=0.2 & Om=uzphi/rsinth+Om0
levm=max([abs(1.-min(Om/Om0)),abs(max(Om/Om0)-1)])
lev=1.-levm+findgen(nlev)/float(nlev-1)*2.*levm
contour,Om/Om0,x,y,/iso,nl=nlev,/fi,lev=lev
colorbar,range=[min(lev),max(lev)],pos=[0.82,0.7,0.84,0.85],/vert,ytickformat='(F5.2)',yticks=2,ytickv=[min(lev),1.,max(lev)],yaxis=0,char=2.5
;
uuh=fltarr(nr,ntheta,3) & uuh(*,*,0)=uxphi & uuh(*,*,1)=uyphi & uuh(*,*,2)=uzphi/x
for ix=0,nx-1 do begin
  for iy=0,ny-1 do begin
    ux[ix,iy]=uuh[ix,iy,0]*sin(theta[iy])+uuh[ix,iy,1]*cos(theta[iy])
    uy[ix,iy]=uuh[ix,iy,0]*cos(theta[iy])-uuh[ix,iy,1]*sin(theta[iy])
  end
end
;
partvelvec,ux,uy,xcoord,ycoord,fraction=frac,/over,length=0.05,col=255
;
!p.charsize=1.4 & xyouts,0.6,0.83,'!6Run A3' & !p.charsize=2.7
;
restore,file='../cool128x256x128a4/stresses.sav'
Om0=0.5 & Om=uzphi/rsinth+Om0
levm=max([abs(1.-min(Om/Om0)),abs(max(Om/Om0)-1)])
lev=1.-levm+findgen(nlev)/float(nlev-1)*2.*levm
contour,Om/Om0,x,y,/iso,nl=nlev,/fi,lev=lev
colorbar,range=[min(lev),max(lev)],pos=[0.15,0.2,0.17,0.35],/vert,ytickformat='(F5.2)',yticks=2,ytickv=[min(lev),1.,max(lev)],yaxis=0,char=2.5
;
uuh=fltarr(nr,ntheta,3) & uuh(*,*,0)=uxphi & uuh(*,*,1)=uyphi & uuh(*,*,2)=uzphi/x
for ix=0,nx-1 do begin
  for iy=0,ny-1 do begin
    ux[ix,iy]=uuh[ix,iy,0]*sin(theta[iy])+uuh[ix,iy,1]*cos(theta[iy])
    uy[ix,iy]=uuh[ix,iy,0]*cos(theta[iy])-uuh[ix,iy,1]*sin(theta[iy])
  end
end
;
partvelvec,ux,uy,xcoord,ycoord,fraction=frac,/over,length=0.05,col=255
;
!p.charsize=1.4 & xyouts,0.6,0.83,'!6Run A4' & !p.charsize=2.7
;
;restore,file='../cool128x256x128a5/stresses.sav'
restore,file='../cool128x256x128a55/stresses.sav'
Om0=1. & Om=uzphi/rsinth+Om0
levm=max([abs(1.-min(Om/Om0)),abs(max(Om/Om0)-1)])
lev=1.-levm+findgen(nlev)/float(nlev-1)*2.*levm
contour,Om/Om0,x,y,/iso,nl=nlev,/fi,lev=lev
colorbar,range=[min(lev),max(lev)],pos=[0.485,0.2,0.505,0.35],/vert,ytickformat='(F5.2)',yticks=2,ytickv=[min(lev),1.,max(lev)],yaxis=0,char=2.5
;
uuh=fltarr(nr,ntheta,3) & uuh(*,*,0)=uxphi & uuh(*,*,1)=uyphi & uuh(*,*,2)=uzphi/x
for ix=0,nx-1 do begin
  for iy=0,ny-1 do begin
    ux[ix,iy]=uuh[ix,iy,0]*sin(theta[iy])+uuh[ix,iy,1]*cos(theta[iy])
    uy[ix,iy]=uuh[ix,iy,0]*cos(theta[iy])-uuh[ix,iy,1]*sin(theta[iy])
  end
end
;
partvelvec,ux,uy,xcoord,ycoord,fraction=frac,/over,length=0.05,col=255
;
!p.charsize=1.4 & xyouts,0.6,0.83,'!6Run A5' & !p.charsize=2.7
;
restore,file='../cool128x256x128a6/stresses.sav'
Om0=2. & Om=uzphi/rsinth+Om0
levm=max([abs(1.-min(Om/Om0)),abs(max(Om/Om0)-1)])
lev=1.-levm+findgen(nlev)/float(nlev-1)*2.*levm
contour,Om/Om0,x,y,/iso,nl=nlev,/fi,lev=lev
colorbar,range=[min(lev),max(lev)],pos=[0.82,0.2,0.84,0.35],/vert,ytickformat='(F5.2)',yticks=2,ytickv=[min(lev),1.,max(lev)],yaxis=0,char=2.5
;
uuh=fltarr(nr,ntheta,3) & uuh(*,*,0)=uxphi & uuh(*,*,1)=uyphi & uuh(*,*,2)=uzphi/x
for ix=0,nx-1 do begin
  for iy=0,ny-1 do begin
    ux[ix,iy]=uuh[ix,iy,0]*sin(theta[iy])+uuh[ix,iy,1]*cos(theta[iy])
    uy[ix,iy]=uuh[ix,iy,0]*cos(theta[iy])-uuh[ix,iy,1]*sin(theta[iy])
  end
end
;
partvelvec,ux,uy,xcoord,ycoord,fraction=frac,/over,length=0.05,col=255
;
!p.charsize=1.4 & xyouts,0.6,0.83,'!6Run A6' & !p.charsize=2.7
;
!p.multi=0
;
;device,/close_file
;set_plot,'X'
;
!p.charthick=1 & !p.thick=1 & !x.thick=1 & !y.thick=1
;
end
